rm(list=ls())

# Set working directory 
setwd("")
source("procedures.R")
library(sandwich)

# Settings
n.large <- 500
n.small <- 100
n.mc <- 10000
n.rw <- 120

# Choose DGP
dgp.nr <- 4

# Random seed
set.seed(22)

# Initialize results matrix
dm.small <- dm.large <- sc.small <- sc.large <- matrix(0, n.mc, 5)

# Parameter estimates
pars <- replicate(n.mc, matrix(0, n.large, 6), simplify = F)

# True and predicted probabilities
prob.pred <- replicate(n.mc, matrix(0, n.large, 6), simplify = F)
prob.true <- replicate(n.mc, rep(0, n.large), simplify = F)

# Load objective functions
of <- getof()

# Loop over MC iterations
for (mm in 1:n.mc){
  
  # Simulate data
  dat <- sim(dgp = dgp.nr, n = n.large + n.rw)
  aux <- matrix(0, n.large, 5)
  
  if (mm %% 1000 == 0) cat(paste("Now at iteration", mm))
  
  # Loop over rolling windows
  for (l in n.rw:(n.large + n.rw - 1)){
    
    # Indexes of estimation sample
    si <- (l - n.rw + 1):l
    
    # Estimate parameter under various scoring rules
    thetas <- rep(0, 6)
    for (j in 1:6){
      thetas[j] <- optimize(of[[j]], y = dat$y[si], x = as.matrix(dat$x[si]), interval = c(-2, 2))$minimum
      pars[[mm]][(l-n.rw+1), j] <- thetas[j]
    }
    
    # Compute loss of MLE minus loss of matching
    for (sr in 2:6){
      aux[(l-n.rw+1), sr-1] <- of[[sr]](b = as.matrix(thetas[1]), y = dat$y[l+1], x = dat$x[l+1], sel = 2) - of[[sr]](b = as.matrix(thetas[sr]), 
                                                                                                                          y = dat$y[l+1], x = dat$x[l+1], sel = 2)
    }	
    
    # True and predicted probabilities
    prob.pred[[mm]][(l-n.rw+1), ] <- logit(thetas*dat$x[l+1])
    prob.true[[mm]][(l-n.rw+1)] <- dat$p.fct(dat$x[l+1])
    
  }
  
  # Compute t-stats for test of EPA
  dm.small[mm, ] <- apply(aux[1:n.small,], 2, dmtest)
  dm.large[mm, ] <- apply(aux, 2, dmtest)
  sc.small[mm, ] <- colMeans(aux[1:n.small,])
  sc.large[mm, ] <- colMeans(aux)
  
  # Save results
  if (mm %% 1000 == 0) save.image(file = paste0("rwMC_dgp", dgp.nr, "_", n.rw, "_", mm, ".RData"))
  
}
